Автор: Винник Екатерина Петровна, 22.М04
import numpy as np
def get_fft(img, shift=False):
h = img.shape[0]
w = img.shape[1]
img = img.astype(complex)
x_arr = np.arange(h)
y_arr = np.arange(w)
x_mat = np.repeat(x_arr.reshape((h, 1)), w, 1)
y_mat = np.repeat(y_arr.reshape((1, w)), h, 0)
if shift:
exp_val = (x_mat + y_mat) * np.pi * 1j
img = np.multiply(img, np.exp(exp_val))
if len(img.shape) == 3:
img = np.swapaxes(img, 0, 2)
img_fft = np.fft.fft2(img)
else:
img_fft = np.fft.fft2(img)
return img_fft
import numpy as np
def get_ifft(img_fft, shift=False):
h = img_fft.shape[0]
w = img_fft.shape[1]
img_fft = img_fft.astype(complex)
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
if shift:
exp_val = (u_mat + v_mat) * np.pi * 1j
img_fft = np.multiply(img_fft, np.exp(exp_val))
if len(img_fft.shape) == 3:
img_fft = np.swapaxes(img_fft, 0, 2)
img_ifft = np.fft.ifft2(img_fft)
else:
img_ifft = np.fft.ifft2(img_fft)
return img_ifft
import numpy as np
def apply_laplac_HP(img_fft):
h = img_fft.shape[0]
w = img_fft.shape[1]
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
dist_mat = np.sqrt(np.power((u_mat - h//2), 2) + np.power((v_mat - w//2), 2))
filter_mat = -4 * np.pi * np.pi * np.power(dist_mat, 2)
if len(img_fft.shape) == 3:
channels = img_fft.shape[-1]
for ch in range(channels):
img_fft[:, :, ch] = np.multiply(img_fft[:, :, ch], filter_mat)
else:
img_fft = np.multiply(img_fft, filter_mat)
return img_fft, filter_mat
import numpy as np
def apply_gauss_HP(img_fft, d0perc=.5):
h = img_fft.shape[0]
w = img_fft.shape[1]
d0 = d0perc * np.min([h, w]) / 2
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
dist_mat = np.sqrt(np.power((u_mat - h//2), 2) + np.power((v_mat - w//2), 2))
filter_mat = np.full_like(dist_mat, 1.) - np.exp(-np.power(dist_mat, 2) / (2 * d0 * d0))
if len(img_fft.shape) == 3:
channels = img_fft.shape[-1]
for ch in range(channels):
img_fft[:, :, ch] = np.multiply(img_fft[:, :, ch], filter_mat)
else:
img_fft = np.multiply(img_fft, filter_mat)
return img_fft, filter_mat
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
img_1 = Image.open('./text1.jpg')
img_2 = Image.open('./text2.jpg')
# represent the images as tensors
np_1 = np.array(img_1)
np_2 = np.array(img_2)
print('Picture dimensions:')
print(np_1.shape)
print(np_2.shape)
Picture dimensions: (256, 256, 3) (256, 256, 3)
# simplify the images by summing the channels
np_1_sum = np.sum(np_1, axis=2) // 3
np_2_sum = np.sum(np_2, axis=2) // 3
np_1_fft = get_fft(np_1_sum, True)
np_2_fft = get_fft(np_2_sum, True)
np_1_gs1, np_1_gs1_fm = apply_gauss_HP(np_1_fft, .01)
np_1_gs2, np_1_gs2_fm = apply_gauss_HP(np_1_fft, .05)
np_1_gs3, np_1_gs3_fm = apply_gauss_HP(np_1_fft, .10)
np_1_gs1i = np.abs(get_ifft(np_1_gs1))
np_1_gs2i = np.abs(get_ifft(np_1_gs2))
np_1_gs3i = np.abs(get_ifft(np_1_gs3))
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].imshow(np_1_gs1i, cmap='gray')
axs[0, 1].set_title('Gaussian highpass filter, D_0 = 1%H')
axs[1, 0].imshow(np_1_gs2i, cmap='gray')
axs[1, 0].set_title('Gaussian highpass filter, D_0 = 5%H')
axs[1, 1].imshow(np_1_gs3i, cmap='gray')
axs[1, 1].set_title('Gaussian highpass filter, D_0 = 10%H')
plt.show()
np_1_ll, _ = apply_laplac_HP(np_1_fft)
np_1_l1i = np.abs(get_ifft(np_1_ll))
fig, axs = plt.subplots(1, 2, figsize=(15, 15), dpi=120)
axs[0].imshow(np_1_sum, cmap='gray')
axs[0].set_title('Original of the 1st image')
axs[1].imshow(np_1_gs1i, cmap='gray')
axs[1].set_title('Laplacian highpass filter')
plt.show()
В результате применения фильтра высоких частот, использующего Лапласиан $D_0 \sim 5\%$ повысилась резкость изображения, текст стал более читаемым. Таким образом, лапласиан в частотной области позволяет добиться того же эффекта, что и лапласиан в пространственной области.
import numpy as np
def apply_amp_HP(img_fft, filter_mat, k1=1., k2=1.):
img_fft = np.multiply((k1 + k2*filter_mat), img_fft)
return img_fft
np_1_hpgs1 = apply_amp_HP(np_1_fft, np_1_gs1_fm, 2, 3)
np_1_hpgs2 = apply_amp_HP(np_1_fft, np_1_gs2_fm, 2, 3)
np_1_hpgs3 = apply_amp_HP(np_1_fft, np_1_gs3_fm, 2, 3)
np_1_hpgs1i = np.abs(get_ifft(np_1_hpgs1))
np_1_hpgs2i = np.abs(get_ifft(np_1_hpgs2))
np_1_hpgs3i = np.abs(get_ifft(np_1_hpgs3))
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].imshow(np_1_hpgs1i, cmap='gray')
axs[0, 1].set_title('Amplification highpass filter (Gaussian, D_0 = 1%H)')
axs[1, 0].imshow(np_1_hpgs2i, cmap='gray')
axs[1, 0].set_title('Amplification highpass filter (Gaussian, D_0 = 5%H)')
axs[1, 1].imshow(np_1_hpgs3i, cmap='gray')
axs[1, 1].set_title('Amplification highpass filter (Gaussian, D_0 = 10%H)')
plt.show()
Резкость изображения в результате усиления высоких частот повысилась, стало возможно прочесть слово susceptibility. На изображении появился незначительный шум.
Изображение может быть представлено в виде произведения его освещенности $i(x, y)$ и коэффициента его отражения $r(x, y)$: $$ \begin{equation} f(x, y) = i(x, y)r(x, y). \tag{2} \end{equation} $$ Также можно представить это произведение в сумму, взяв логарифм от функции $f(x,y)$ : $$ \begin{equation} z(x, y) = \ln f(x, y) = \ln i(x, y) + \ln r(x, y), \tag{3} \end{equation} $$, где $Z(u, v) = \mathfrak{F}[z(x, y)]$ сумма $\ln i$ и $\ln r$ в частотной области. Обозначим $\mathfrak{F}[\ln i(x, y)]$ как $F_i(u, v)$. $F_r(u, v)$ аналогично соответствует $r(x, y)$. Отсюда, $$ \begin{equation} Z(u, v) = F_i(x, y) + F_r(x, y). \tag{4} \end{equation} $$ В результате применения фильтра $H(u, v)$, можем получить преобразованное изображение в пространственной области следующим образом: $$ \begin{equation} g(x, y) = e^{s(x, y)} = e^{i'(x, y)} e^{r'(x, y)}, \tag{5} \end{equation} $$ где $i'(x, y) = \mathfrak{F}^{-1}[H(u, v)F_i(u, v)]$ and $r'(x, y) = \mathfrak{F}^{-1}[H(u, v)F_r(u, v)]$.
Будем использовать модифицированный фильтр Гаусса как гомоморфный: $$ \begin{equation} H(u, v) = (\gamma_H - \gamma_L)\left[ 1 - e^{-c(D^2(u, v)/D_0^2)} \right] + \gamma_L, \tag{6} \end{equation} $$ где $c$ и $D_0$ константы, $D(u, v)$ есть расстояние от точки $(u, v)$ до центра изображения, $\gamma_H > 1$ and $\gamma_L < 1$ -- параметры, регулирующие то, как фильтр воздействует на высокочастотные (соответствующие отражениям) и низкочастотные (соответствующие внешнему освещению) части изображения.
import numpy as np
def apply_gs_hom(img_fft, g_h, g_l, c=1., d0perc=0.5):
h = img_fft.shape[0]
w = img_fft.shape[1]
d0 = d0perc * np.min([h, w]) / 2
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
dist_mat = np.sqrt(np.power((u_mat - h//2), 2) + np.power((v_mat - w//2), 2))
filter_mat = (g_h - g_l) * (np.full_like(dist_mat, 1.) - np.exp(-c*np.power(dist_mat, 2)/(d0*d0))) + \
g_l * np.full_like(dist_mat, 1.)
img_fft = np.multiply(img_fft, filter_mat)
return img_fft
# make z(x, y)
z_1 = np.log((np_1_sum / 255))
z_2 = np.log((np_2_sum / 255))
Применим БПФ к изображениям:
z_1_fft = get_fft(z_1, True)
z_2_fft = get_fft(z_2, True)
# log of abs
z_1_fft_logabs = np.log2(np.abs(z_1_fft)+1)
z_2_fft_logabs = np.log2(np.abs(z_2_fft)+1)
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].matshow(z_1_fft_logabs, cmap='gray')
axs[0, 1].set_title('Log of abs. of FFT of the 1st image in the z-form')
axs[1, 0].imshow(np_2_sum, cmap='gray')
axs[1, 0].set_title('Original of the 2nd image')
axs[1, 1].matshow(z_2_fft_logabs, cmap='gray')
axs[1, 1].set_title('Log of abs. of FFT of the 2nd image in the z-form')
plt.show()
Применим гомоморфный фильтр для различных $\gamma_H$ и $\gamma_L$ к первому изображению.
# D0 ~ 5%, different gamma and c
s_fft_gshm_1_1 = apply_gs_hom(z_1_fft, 1.5, 0.5, 1, 0.05)
s_fft_gshm_1_2 = apply_gs_hom(z_1_fft, 1.5, 0.5, 5, 0.05)
s_fft_gshm_1_3 = apply_gs_hom(z_1_fft, 2, 0.25, 1, 0.05)
s_fft_gshm_1_4 = apply_gs_hom(z_1_fft, 2, 0.25, 5, 0.05)
# D0 ~ 50%, different gamma and c
s_fft_gshm_2_1 = apply_gs_hom(z_1_fft, 1.5, 0.5, 1, 0.5)
s_fft_gshm_2_2 = apply_gs_hom(z_1_fft, 1.5, 0.5, 5, 0.5)
s_fft_gshm_2_3 = apply_gs_hom(z_1_fft, 2, 0.25, 1, 0.5)
s_fft_gshm_2_4 = apply_gs_hom(z_1_fft, 2, 0.25, 5, 0.5)
s_1_1 = np.abs(get_ifft(s_fft_gshm_1_1))
s_1_2 = np.abs(get_ifft(s_fft_gshm_1_2))
s_1_3 = np.abs(get_ifft(s_fft_gshm_1_3))
s_1_4 = np.abs(get_ifft(s_fft_gshm_1_4))
s_2_1 = np.abs(get_ifft(s_fft_gshm_2_1))
s_2_2 = np.abs(get_ifft(s_fft_gshm_2_2))
s_2_3 = np.abs(get_ifft(s_fft_gshm_2_3))
s_2_4 = np.abs(get_ifft(s_fft_gshm_2_4))
g_1_1 = np.exp(s_1_1)
g_1_2 = np.exp(s_1_2)
g_1_3 = np.exp(s_1_3)
g_1_4 = np.exp(s_1_4)
g_2_1 = np.exp(s_2_1)
g_2_2 = np.exp(s_2_2)
g_2_3 = np.exp(s_2_3)
g_2_4 = np.exp(s_2_4)
fig, axs = plt.subplots(4, 2, figsize=(15, 30), dpi=120)
axs[0, 0].matshow(g_1_1, cmap='gray')
axs[0, 0].set_title('D0 ~ 5%, g_h = 1.5, g_l = 0.5, c = 1')
axs[0, 1].matshow(g_1_2, cmap='gray')
axs[0, 1].set_title('D0 ~ 5%, g_h = 1.5, g_l = 0.5, c = 5')
axs[1, 0].matshow(g_1_3, cmap='gray')
axs[1, 0].set_title('D0 ~ 5%, g_h = 2, g_l = 0.25, c = 1')
axs[1, 1].matshow(g_1_4, cmap='gray')
axs[1, 1].set_title('D0 ~ 5%, g_h = 2, g_l = 0.25, c = 5')
axs[2, 0].matshow(g_2_1, cmap='gray')
axs[2, 0].set_title('D0 ~ 50%, g_h = 1.5, g_l = 0.5, c = 1')
axs[2, 1].matshow(g_2_2, cmap='gray')
axs[2, 1].set_title('D0 ~ 50%, g_h = 1.5, g_l = 0.5, c = 5')
axs[3, 0].matshow(g_2_3, cmap='gray')
axs[3, 0].set_title('D0 ~ 50%, g_h = 2, g_l = 0.25, c = 1')
axs[3, 1].matshow(g_2_4, cmap='gray')
axs[3, 1].set_title('D0 ~ 50%, g_h = 2, g_l = 0.25, c = 5')
plt.show()
Лучший результат достигается при значениях $D_0 \sim 50\%$, $\gamma_H = 2$, $\gamma_L = 0.25$, $c = 1$. Применим негатив и степенное преобразования, чтобы осветлить полученное в результате гомоморфной фильтрации изображение.
g_2_3_negpow = np.power((1 - g_2_3/np.max(g_2_3)), 4)
fig, axs = plt.subplots(1, 2, figsize=(12, 6), dpi=120)
axs[0].imshow(np_1_sum, cmap='gray')
axs[0].set_title('Original of the 1st image')
axs[1].matshow(g_2_3_negpow, cmap='gray')
axs[1].set_title('D0 ~ 50%, g_h = 2, g_l = 0.25, c = 5 \n after the negative and power transformations')
plt.show()
В результате применения гомоморфной фильтрации получили наиболее читаемое изображение по сравнению с изображениями, полученными ранее.
# D0 ~ 50%, different gamma and c
s_fft_gshm_3_1 = apply_gs_hom(z_2_fft, 1.5, 0.5, 1, 0.5)
s_fft_gshm_3_2 = apply_gs_hom(z_2_fft, 1.5, 0.5, 5, 0.5)
s_fft_gshm_3_3 = apply_gs_hom(z_2_fft, 2, 0.25, 1, 0.5)
s_fft_gshm_3_4 = apply_gs_hom(z_2_fft, 2, 0.25, 5, 0.5)
s_3_1 = np.abs(get_ifft(s_fft_gshm_3_1))
s_3_2 = np.abs(get_ifft(s_fft_gshm_3_2))
s_3_3 = np.abs(get_ifft(s_fft_gshm_3_3))
s_3_4 = np.abs(get_ifft(s_fft_gshm_3_4))
g_3_1 = np.exp(s_3_1)
g_3_2 = np.exp(s_3_2)
g_3_3 = np.exp(s_3_3)
g_3_4 = np.exp(s_3_4)
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].matshow(g_3_1, cmap='gray')
axs[0, 0].set_title('D0 ~ 50%, g_h = 1.5, g_l = 0.5, c = 1')
axs[0, 1].matshow(g_3_2, cmap='gray')
axs[0, 1].set_title('D0 ~ 50%, g_h = 1.5, g_l = 0.5, c = 5')
axs[1, 0].matshow(g_3_3, cmap='gray')
axs[1, 0].set_title('D0 ~ 50%, g_h = 2, g_l = 0.25, c = 1')
axs[1, 1].matshow(g_3_4, cmap='gray')
axs[1, 1].set_title('D0 ~ 50%, g_h = 2, g_l = 0.25, c = 5')
plt.show()
g_3_3_negpow = np.power((1 - g_3_3/np.max(g_3_3)), 2)
fig, axs = plt.subplots(1, 2, figsize=(12, 6), dpi=120)
axs[0].imshow(np_2_sum, cmap='gray')
axs[0].set_title('Original of the 2nd image')
axs[1].matshow(g_3_3_negpow, cmap='gray')
axs[1].set_title('D0 ~ 50%, g_h = 2, g_l = 0.25, c = 5 \n after the negative and power transformations')
plt.show()
В результате применения гомоморфной фильтрации получили наиболее читаемое изображение по сравнению с изображениями, полученными ранее.